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Abstract 

This paper examines uniqueness and stability results for an inverse problem in thermal 
imaging. The goal is to identify an unknown boundary of an object by applying a heat flux 
and measuring the induced temperature on the boundary of the sample. The problem is 
studied both in the case in which one has data at every point on the boundary of the region 
and the case in which only finitely many measurements are available. An inversion procedure 
is developed and used to study the stability of the inverse problem for various experimental 
configurat ions . 
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1 Introduction 


Thermal imaging is a technique of wide utility in non-destructive testing and evaluation. The 
technique is used to recover information about the internal condition of an object by applying 
a heat flux to its boundary and observing the resulting temperature response on the object’s 
surface. From this information one attempts to determine the internal thermal properties of 
the object, of the shape of some unknown portion of the boundary. 'Thermal imaging has 
been much investigated as a method for detecting damage or corrosion in aircraft. See [10] 
for an account of the technology and typical data processing techniques that are employed, 
and a more extensive bibliography on the subject. 

One of the most common uses of thermal imaging is for the detection so-called “back 
surface” corrosion and damage. Briefly, one attempts to determine whether some inaccessible 
poition of an object s boundary has corroded, and therefore changed shape. In this paper 
we investigate the inverse problem of determining changes in the boundary profile of a two- 
dimensional sample by using thermal imaging. We consider a certain portion of the surface 
of a rectangular sample to be accessible for measurements and the remainder of the surface, 
which may be corroded, inaccessible. This problem has been considered by others [3, 4] with 
an emphasis on recovering estimates of the unknown surface from data by using an output 
least-squares method. 

We examine both a continuous and finite data version of the inverse problem. The 
continuous version assumes that one has data at every point on the accessible portion of the 
object, s surface. The finite data version assumes that only finitely many measurements have 
been made. Our goals are 

• To examine uniqueness and continuous dependence results for the continuous version 
of the inverse problem, and what they imply for the finite data inverse problem. 

• To examine how various experimental parameters affect stability and resolution for the 
finite data inverse problem, especially the effect of measurement locations on stability. 

• To determine how one might incorporate a priori information or assumptions into the 
finite data inverse problem. 
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Our main focus is not to develop inversion algorithms, but in the course of examining the 
problem we derive (but do not prove convergence for) an inversion procedure for the finite 
data inverse problem. This algorithm allows the easy incorporation of a priori assumptions 
into the inversion process. We apply the algorithm to several simulated data sets to illustrate 
our conclusions. Our study of the stability of the inverse problem reduces to studying the 
invertibility of a certain matrix, which we do with a singular value decomposition. We do 
not make any explicit finite dimensional parameterization of the unknown surface. 

We should note that a very similar approach has been used in [8] to study resolution and 
stability for the inverse conductivity problem. Isaacson, Cheney, and others [6, 7] have also 
carried out similar sensitivity studies related to the inverse conductivity problem, especially 
the effect of finite^ many measurements on the inversion process. 

The outline of the paper is as follows. In Section 2 we present the mathematical for- 
mulation of the continuous and finite data versions of the inverse problem. In Section 3 we 
derive a linearized version of the continuous problem, and in Section 4 we show how this 
leads (as thermal inverse problems often do) to a first kind integral equation which must be 
inverted. In Section 5 we use the integral equation formulation to examine uniqueness and 
stability lesults for the linearized version of the inverse problem. In Section 6 we consider an 
algoiithm for solving the finite data version of the inverse problem and how this approach 
can be used quantify the stability of the problem. In the last section we present a variety 
of numerical studies to examine the effects that various experimental parameters have on 
the stability and resolution of the inversion process, and the effect of incorporating a priori 
assumptions into the inversion procedure. 


2 The Inverse Problem 

Consider a sample to be imaged as a two-dimensional region Q lying between the two surfaces 
x 2 = S(x i) and x 2 = 1 as illustrated below. 
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x 2 = 1 


I x 2 = S(x, ) 

X 

1 

Figure 1: Sample geometry. 

We will refer to the surface x 2 = 1 as the “top” or “front” surface and x 2 — S(x i) as the 
“back” surface. We assume that the ends of the sample are sufficiently far away that they 
can be ignored, so for our purposes the sample is unbounded in the x\ direction. The top 
surface is assumed to be accessible for inspection or measurements, but the back surface 
x 2 = S{x i) is inaccessible. This is the portion of the sample to be inspected for corrosion. 
The ideal uncorroded case is a flat, back surface S(x 1 ) = 0. In the corroded case above 
S(xi) > 0 for some values of x±. We will assume that the function S belongs to iZ 2 (IR), 
although this assumption will later be relaxed. In particular, since H 2 (JR) c C^IR) there is 
a continuous unit normal vector field on the back surface. The goal is to determine the back 
surface or the function S by taking measurements only on the front surface. 

A time-dependent heat flux g(x\,t) is applied to the top of the sample x 2 — 1. We will 
assume that the sample material is homogeneous with thermal diffusivity k and thermal con- 
ductivity a, both known constants. We will use T(x , t) to denote the resulting temperature 
induced in Q, where x = (x 1 ,x 2 ). The direct thermal diffusion problem will be modeled by 
the standard heat equation 

= 0 in SI, (2.1) 

= g(xi,t) on x 2 = 1, 

— 0 on x 2 = S(x i), 

= T 0 (x), 

for t > 0, where denotes the outward normal derivative on the boundary of SI. The 
function T 0 (x) is the initial temperature of the region O at time t = 0. Note that the back 
surface is assumed to block all heat conduction. 


dT 

at 


-k AT 


a 


dT 

dv 

dT 


a - 


dv 
T(x, 0) 
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We will also assume that the heat flux g(x 1 ,t) is periodic, of the form Re[#(.ri)e* wf ] with 
uj > 0. For simplicity, we also take the constants k and a equal to one. Under these 
assumptions the solution to equation (2.1) is given as T(xJ) = Re[e^u(.r)] where u(x) 
satisfies 

Aa — icon = 0 in (2.2) 

dn 

^ = #04 ) on x 2 = 1, 
du 

— = 0 on x 2 = S(xi), 

at least after transients from the initial condition have died out. The main case of interest is 
that in which g(x i) is constant, corresponding to uniform heating of the outer surface. This 
is typically the case when heat or flash lamps are used to provide the input flux g. For the 
moment, however, we will not restrict g. 

There are two versions of the inverse problem to be considered: 


Continuous Version: 

Given measurements of u(x) at all points on the top surface x 2 = 1, determine the 
function 5'(;r 1 ). 

Finite Data Version: 

Given measurements of u{x) on the top surface x 2 = 1 at points x x = a 1} a 2 , . . . , a„ , 
estimate the function i>(;ri). 

The finite data version corresponds to the case in which one has actual measurements. The 
data need not be actual point measurements of the temperature «, but this is the most 
common situation. Of particular interest are the questions 

1. Can the function S(x-i) be uniquely determined by knowing u(x{) for all x± on the top 
surface? 

2. If S'(.ri) is uniquely determined by u(aq), how sensitive is ^(.rj) to perturbations in the 
data? What kinds of features in the back surface x 2 = ^(oq) can or cannot be easily 
determined from the data? 
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3. Since any practical application falls under the finite data formulation, how stable is 
the estimate of S(xi) based on finitely many pieces of data? What factors influence 
stability in this case, and is there an inversion procedure to produce a reasonable 
estimate of S(x i) using finitely many measurements? 

The first question is easily answered “yes” by a standard argument. This is the content of 
the following -result. . - 

Theorem 2.1 (Uniqueness) Let u(xi,x 2 ‘, S) denote the solution to (2.2) with back surface 
S and nonzero flux g. If u(xi, 1; Si) = «(;i'i,l;S 2 ) for each (;ri,l) in an open subset C of 
the top surface of f 2 , then Si = S 2 . 

Proof: Suppose S 2 ^ S 2 . Using the shorthand notation m = u(x 1 ,x 2 ;S i ), we have that u x 
and u 2 have the same Cauchy data on C, and by unique continuation ui and u 2 agree on any 
connected domain on which both are defined, provided that that domain contains an open 
portion of C . Assume that Si and S 2 are not equal and denote by the region bounded 
bv x 2 = 1 and x 2 = Sflx i). Si ^ S 2 , so there is some non-empty connected component 
Z> of Hi \ 0 2 or n 2 \ Ui. Let us suppose the latter so that the region D is bounded by 
X ‘2 = Si (xi) and x 2 = S 2 (xi). On x 2 = S 2 (x i) we know that the normal derivative of u 2 is 
identically zero; on x 2 = Si(aq) we know that the normal derivative of u\ is zero (from inside 
fU) and since u 2 = Ui and u 2 is smooth across x 2 = Si(.ri), we conclude that the normal 
derivative of u 2 vanishes on the boundary of D. This forces u 2 = 0 inside D. Standard elliptic 
regularity arguments force u 2 = 0 inside ft 2 , which implies that the flux g is identically zero, 
a contradiction. Thus we must have Si = S 2 and so the back surface S (.Ti ) is uniquely 
determined by the boundary data on am' open portion of the top surface. ■ 

The second and third questions will be examined in the next few sections by considering 
a linearization of the original problem. 


3 A Linearization 

We now linearize the original direct problem given by equation (2.2) with respect to the 
function S. Let Wo(>^ : ) denote the solution to (2.2) with S = 0. The surface x 2 ~ 0 is the 
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point at which we will linearize, since this represents the uncorroded or ideal profile from 
which we hope to detect any deviation. 


Let O 0 denote the region {(.Ti,.t 2 ) : 0 < x 2 < 1} in 1R 2 . For any 5 G Ff 2 (R) C C^IR), 
let us construct the map <p from Q to fl 0 by 


d>(x\,x 2 ) = (:r : , 


,r 2 - 5(>ri) 

1 - S( Xl ) } 


If |5(xi)| < 1 then it is easy to check that cp is invertible on fi 0 - 


This map fixes the top 


boundary of f2 and maps the bottom surface to x 2 = 0. Let y = <b(x) and v(y ) = u(<p~ l (y)) = 
u(x). Under such a change of coordinates V x = ( D(p) T V y and JL- = JL where r/ y = ( Dcp)v x , 
(D<p) is the derivative of 0, and (D<p) T is the transpose. Under this change of coordinates 
the boundary value problem (2.2) becomes 


V • kVv — iuv = 0 in fi 0 , (3.3) 

dv 

foj = ^(^i) on x 2 = 1, 
dv 

— = 0 on X‘ 2 = 0, 

07 ] 

where k = (Dd>)(D(p) T . The vector y can be written as 


0 

, i 

on To — 1 n — 

S'(afi) 

1 

V i ’ f j • 

VI + (5'(.n)) 2 

l + (S'(x ,)) 2 

L 5(an)-l J 

L 1-5(X!) J 


This change of coordinates removes the unknown 5 from the definition of the boundary of 
O and puts S into the coefficients of the heat operator and boundary conditions. 

Now we linearize the problem with respect to S about S = 0 by assuming that 5 = eS 
for some function 5, where e is some small real number. Let u denote the solution to (3.3) 
for a general S and iiq the solution to equation (3.3) for ,5 = 0. Suppose that the resulting 
perturbation in u G can be represented as u = u 0 + €V o- If we substitute these relations into 
(3.3), use the fact that uq satisfies (3.3) with S — 0, and drop all quadratic e terms then we 
obtain a linearized version of the direct problem, 


At’o — itov o 

dvo 

dv 
dv o 
dv 


— V-[/cVm 0 ], in n 0 , 
-S{yi)9(yi), on y 2 = 1 
—S'{yi)d Q {y\), on y 2 = 0, 


(3.4) 
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where u,(y) is the matrix 


&(y) = 


0 (y 2 - l)S'{yi) 


[ ( 2/2 - 1)5'(2/i) 25(2/0 J 

In particular, we are interested in the special case g = 1 corresponding to spatially 
uniform periodic heating of the top surface. Under this condition the function Uq depends 
only on x 2 and equation (3.3) for u 0 becomes a two-point boundary value problem 


itQ — iunio — 0 on (0,1), (3.5) 

«i(0) = 0, 

w o(l) = 1- 

If we now use u rather than i?o to denote the temperature perturbation satisfying equation 
(3.3) then the linearized version of the problem for g = 1 is then 


Au-ium = (1 -x 2 )u' 0 (x 2 )S r, (xi) - 2u'£(x 1 )S(x 1 ) on O 0 , (3.6) 

du 

— = —S(xi) on x 2 = 1, 

du 

TT~ = 0 Oil X 2 = 0. 

du 

For simplicity, this is the version of the problem we will examine, although the more 
general linearized version (3.4) can be examined using similar techniques. Note that the 
linearized problem is defined on the domain O 0 which does not depend on S. 


4 An Integral Identity 

Let the function d(a ) = u(a, 1) denote the top surface data from the direct problem (3.6). 
Given that the relation between S and d is linear, it seems reasonable that this relationship 
can be expressed by an integral operator 

/ OO 

4>a{y)S(y) dy 

-oo 

where 6 a is some function which depends on the particular point a. Such a relationship does 
exist and we can say quite a bit about the function(s) <b a , as we now demonstrate. 
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Let L = A — iu). By Green’s second identity 



where d> is any sufficiently regular function defined on fh Assume that u is a solution to 
(3.6) and that <b is a function which satisfies Ld> = 0 on Q 0 and = 0 on the surface x 2 = 0. 
Then the above equation becomes 

[ -X 2 )u' 0 (x 2 )S"(xi) -2ug{x 1 )S'(x 1 ))dx + f°° 4>(x 1 ,l)S(x 1 )dx 1 

J —oo 

deb, 

Let us now complete the specification of <fi by requiring ^ = S a on the top surface, where S a 

denotes a delta function on the surface x 2 — 1 at the point x\ — a. If we write ef) a to denote 

the dependence of ej> on a then we obtain 

f r oo 

4>a(^)((i-x 2 )u' 0 (x 2 )S"(x 1 )-2u'i(x 1 )S(x 1 ))dx+ <p a (xul)S(x 1 )dxi = -d(a). (4.7) 

Note that this equation involves no unknown quantities except S on the left side. 

It is worth saying a few' words about the function eb a w'hich satisfies 


A0 a — ito<f) a — 0 in f2 0 , 

d&a 


(4.8) 


dv 

d4>a 

dv 


= S a on x 2 — 1, 
= 0 on x 2 — 0. 


Let T(x) be a Green s function for the operator A — iev on IR 2 : such a function is given by 


r(.r) = ^(ker(ry/w) + ikei(ry/w )) , 

where r = |.r| = ^Jx\ 4- x 2 and ker() and kei() are the Kelvin functions (see [1], section 9.9). 
The function ker(r) has a — ln(r) singularity as r approaches zero, while kei(r) is bounded. 
Both functions and their derivatives are smooth away from zero and rapidly decreasing as r 
tends to infinity, where “rapidly decreasing” means faster than any power of A If we define 
W a {'?) = ~~2F(x — x a ) where x a is the point (a, 1) on the top surface, then standard potential 
theory arguments ([9], chapter 3) show that 


Alp a — i.UJlpa 
dj’g 
dv 


0 in O 0 , 

S a on x 2 — 1. 


(4-9) 
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It is not true that 


However, if we take v a G 


dii’c 


dv 

H 1 ( f2) to satisfy 


= 0 


on X 2 = 0. 


Av a — i(jjv a 
dvg 
dv 


dv 


— 0 in f2o, 


= 0 on X 2 — 1, 


dtpa. 

dv 


on X 2 = 0. 


(4.10) 


then 4> a — xp a + v a . Since the Neumann data for v a on the bottom surface is — which 
is in H°°{ IR.) (the singularity for ip a lies on the top surface, and away from this singularity 
4' a is smooth and rapidly decreasing), one can show that v a is in H°°(Vl 0 ). As a result, the 
function d> a {x) has a — ~ln|.r| singularity near x — 0 and otherwise is smooth and rapidly 
decreasing in |;r|, along with its derivatives of all orders. 

If we write the integrals in equation (4.7) with limits, we find that S must satisfy 

/ oo 

(p a (,Ti)S"(,Ti) + ga^iqjS^Ti)) dx-i = -d(a) 

-OO 

where 

Pa(^l) = I* <t>a(Xl,X2)u' 0 (X2){l ~ X 2 ) dx 2 , (4.11) 

Qa(x 1 ) = -2 f <b a (x i , x 2 ) u'q (x 2 ) dx 2 + <b a 1 , 1 ) . (4-12) 

J 0 

One can check that the integral in q a (x i) is continuous as a function of x\, smooth away 
from X\ — «, and rapidly decreasing in xi. Also, since 4> a {xi, 1) has a logarithmic singularity 
at = a, so does q a . Moreover q a is an L 2 function. The function p a (x i) is also clearly 
smooth away from ,ri = a and rapidly decreasing in xi. The singularity of 4> a (x i,a' 2 ) looks 
like the singularity of ker|.r — ar a |, and one can use this fact to expand the function <b a (x !, rr 2 ) 
near the singularity to show that the function p a (x i) is actually in 7f 2 (IR). Since both p a (: r x ) 
and q a {x i) tend rapidl}^ to zero as |.Ti| — > oo we can integrate by parts twice to find that 

/ oo />oo 

Pa(#i)S"(.r 1 ) dxi = / p"(xi)S(xi) dxi. 

-oo J — OO 

Equation (4.9) can now be wuitten 


/ OO 

c a (xi)S{xi) dxi - d(a) 

-OO 


(4.13) 
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where 

Ca{zi) = -p"(x 1 ) - q a (x i) <E T 2 (1R)- (4.14) 

Given the translation invariance of this problem in the rri direction and the fact that the 
flux g{xi) = 1 does not depend on X \ , it is clear that c a (.ri) = c{x\ — a) for some function c. 

Lemma 4.1 The function 0 ( 04 ) — M O (O)ic^0o(xi,O), where 4>o satisfies (4-8) with a = 0 and 
u 0 satisfies (3.5). 

In proving this result, We will make use of the following simple fact. 

Remark 4.1 If f(t) and g(t) are functions defined on [0,1] with f{ 0) = /'( 1) = 0 and 
g'( 0) = 0, g'(l) = 1, and g" = iug then 

f f"(t)g'(t)(l ~t)dt = f( 1) - 2iw f 1 f (t)g(t) dt + iwf(0)g(0) + ico 1 - t)f(t)g'(t) dt. 

Jo Jo 

To see this identity requires only a few applications of integration by parts, f udv — uv — 
I vdu. Take u = g'(t)( 1 - t), dv = f'(t) dt to obtain 

£ /"(%'(*)( 1 -t)dt = £ f(t)g'(t) dt - iu £ f'(t)g(t)( 1 - t) dt, 

where we have made use of g" = icvg. Integrate the first integral by parts with u = g' and 
dv ~ f dt: integrate the second integral by parts with u = (1 -t)g and dv = f dt. After a 
number of cancellations, the Remark follows immediately. 

Proof of Lemma 4.1: From equations (4.11), (4.12) and (4.14) we obtain 

c(-ri) = - f () ((! - x 2 )~^u' 0 (x 2 ) - 2(/> 0 Uo(x 1 )) dx 2 - fio(xi, 1 ). 

Recall that u 0 is a function of x 2 only, and that u’( = icjuo with u[,(0) = 0 and m{,( 1 ) = 1 . 
Since A <£> 0 ~ icvfio = 0 we have = itud>o — With this substitution and u { J = iujug 

ri d 2 d> 

c(.ri) = — ^ [iu>( 1 — x 2 )<hou'o(x 2 ) — (1 — ;r 2 )u{ } — ^ 2 ° — 2?A<A)Uo] dx 2 — 4>o{xi, 1). 

To finish the proof, use Remark 4.1 above with f = d> 0 and g = u 0 . After a few cancellations 
the statement of the lemma follows.* 
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Lemma 4.1 implies a number of things about the function c(x i), in particular, that c is 
in H°°(WC). It is important to note that «o(0) is never zero, for uq satisfies the two-point 
boundary value problem (3.5) and one easily finds that 

&<XX2 _|_ g — CtX-2 


a(e c 




Uo(*2) = 

where a = (1 + i) \Juj f2. In particular, 

«o( 0 ) = — 

a(e a — e ~ a ) 

The numerator above is never zero, hence Wq(0) 7 ^ 0. As a result the function c(:ri) will not 


be identically zero for any uj > 0 . 

There is one more fact which will be extremely useful. Defining the Fourier transform 
fly) of a function f(x) of a single variable by 


- f°° 

f(y) = / e~ ixy f(x)dx, 

J — OO 

we have 


Lemma 4.2 For c given by Lemma 4-1, 


c{y) = 


2iojuo(0) 
a{e Q — e~ a ) 


where a = sjy 1 + uj. Also, the function c(y ) is never equal to zero. 


Proof: Let <£ 0 ( 2 / 1 , ^ 2 ) denote the Fourier transform of <£ 0 (^i, x 2 ) with respect, to Xi only, 
where <£ 0 satisfies equation (4.8) with a = 0. Then 

vs r 00 

<£ 0 ( 2 / 1 , ^ 2 ) = / (t>o(x 1 ,x 2 )e~ tXiyi dx x . 

If <£ 0 is sufficiently rapidly decreasing (as in our case) then ^ and the Fourier transform 
operator commute. Fourier transforming both sides of the boundary value problem (4.8) with 
respect to X\ yields a two-point boundary value problem for <£ 0 ( 2 / 1 , ^ 2 ) in the x 2 variable, 

= 0, 

= 1 , 

= 0 , 

(4.15) 


(f 2 4> 0 

dxo 


- (yj + iu)4> 0 


dfh 
dx 2 
df> 0 

dx 2 


(1) 

(0) 
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where we have used 5 0 — 1. The solution is 


4>o(yi,*2) = 


I „—aX2 


e- + e 


a(e a — e~ a ) ' 

where a = \Jvl + iu ;. On the bottom surface x? = 0, 

4>o(yu 0 ) = - 2 

a\e a — e Q ) . - 

which, combined with Lemma 4.1 yields the conclusion. The fact that c(y) is never zero 
follows from lj > 0 and w o (0) ^ 0. a 


5 The Linearized Inverse Problem 


In this section we will examine the continuous version of the linearized inverse problem. 
Suppose that we have top surface data = w(&’i,l) where u satisfies equation (3.6). 

Based on equation (4.13) and the above noted fact that c a (:ti) = c(x 1 — a) we conclude that 
S and d satisfy the relationship 


/ OO 

S{yi)c(yi - aq) dy Y ~ d(aq) 

-OO 


or 

/ OO 

- yi) dyi = S * <£> — d(x{) (5.16) 

-OO 

where 4>(xi) = c( :zq) and denotes convolution. With d{x\) considered known this 
becomes a first kind integral equation for the unknown function S . The kernel <f> is known, 
or can be determined, according to Lemma 4.1. First kind integral equations have been 
extensively studied ([11, 12]) and are well-known to be unstable; small perturbations in the 
right hand side d(x) can lead to arbitrarily large changes in the solution 5. However this 
humiliation of the inverse problem as an integral equation will allow us to obtain uniqueness 
and stability estimates for the linearized version of the inverse problem. 


Theorem 5.1 (Uniqueness) If the data d(x) is an L 2 function and if there exists a solution 
S G L 2 (1R) to the linearized inverse problem, then S is unique. 
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Proof: Suppose Si and S 2 are L 2 functions which both give rise to data d. Let S = Si — S 2 . 
Linearity implies that the data for S is identically zero. Fourier transform both sides of the 
equation S * 4> — 0. By the convolution theorem and the fact that 4> = c if <p(x) = c(—x) 
(where the overbar denotes complex conjugation) we obtain 4>S = 0. By Lemma 4,2 i — 
c 0 and we conclude that 

5 = 0. 


so that 5 = 0 or Si = S 2 . ■ 


Remark 5.1 In the preceding proof, we assumed a priori that S is in L 2 (IR). In general . 
for an arbitrary d E L 2 (3R) we cannot find a function S in L 2 which gives rise to data d via 
equation (5.16). 


The convolution equation (5.16) also provides information on continuous dependence. Since 
the function f>o is smooth and never equal to zero, we can define the space of functions L 2 (IR) 
with the norm 



From Lemma 4.2 it follows that grows like ze*. The norm || ||* thus puts a heavy 

penalty on high frequencies; the functions in this space are very smooth. We can Fourier 
transform both sides of equation (5.16), divide by 0 O and take the L 2 norms of both sides to 
obtain 

Theorem 5.2 (Continuous Dependence) If a back surface x 2 = S(x 1 ) generates front sur- 
face data d(x) for the linearized problem (3.6), then 

IISIU* < c\\d\u 

where C is independent of d. 

Estimates of S from data d will thus be extremely sensitive to any noise, because the in- 
version process weights a frequency / in the data by a factor proportional to fe f . Lemma 4.2 
and the structure of the convolution operator mapping S to the data d make it clear that it 
will be difficult to estimate the high spatial frequency components in the Fourier decompo- 
sition of S , for these components are heavily damped out by the forward mapping. 
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6 The Case of Finitely ]Vf any ]VIeasurements 


Suppose that, we have point estimates rf(a,-) = w(a t - , 1) of the temperature on the top surface 
at n distinct points. How can we construct a reasonable estimate of the function S(xi)? 
How can we quantify the stability of the reconstruction with respect to errors in the data, 
and how' does the choice of measurement locations a z - affect the stability? Let us assume 
that we seek an estimate S € L 2 (1R). Physical considerations make it desirable to obtain an 


estimate with moie regularity, but this will be a consequence of the proposed reconstruction 
procedure. Based on the convolution equation (5.16) we know that S' must satisfy the n 
constraints 


< 


/ OO 

Sfyi)c z (.Ti) dxi =d(di ), 

-CO 


* = 1 ,. 


n, 


(6.17) 


with Qfyi) = c(ai — xi) where cfyq) is the function from Lemma 4.1 and < /, g >= f u fg is 
the usual L 2 inner product. Note that since c* is an L 2 function, S M-< S, c* > is a bounded 
linear functional on L 2 . The set (6.17) is a horribly underdetermined set of equations. We 
can expect to find an entire translated subspace of functions of codimension n in L 2 (IR) 
which satisfy the given conditions, and any such function “solves” the inverse problem, in 
the sense that it gives rise to the measured data. 

One practical method for specifying a unique function in L 2 which solves the inverse 
problem is to seek that element in L 2 which satisfies the given conditions and has minimal 
norm. That such an element exists follows from the fact that the relations (6.17) define a 
closed convex subset of L 2 and hence this subset has a unique element of minimal norm . This 
idea has been used before ([8]) to construct a “pseudo-inverse” for the finite measurement 
case and to characterize the stability and information content for the inverse conductivity 
problem, and has also been used for reconstruction from partial information in tomographic 
problems [5]. The approach has several merits: In the present case it leads to an exceptionally 
easy and efficient inversion algorithm which allows us to study the conditioning of the inverse 
problem independent of any explicit finite dimensional parameterization of the unknown S. 
By weighting the L 2 space appropriately we can also incorporate a priori assumptions into 
the reconstruction procedure and examine the effect these assumptions have on stability. 
Also, given the continuous dependence result from Lemma 5.2 and the fact that data is 
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invariably noisy, we know that any inversion procedure will tend to give extraneous high 
frequency components in an}' estimate of S ; choosing the estimate of minimal norm should 
help to damp out spurious components in the estimate. In this sense the procedure may be 
viewed as a form of regularization. 

It is an easy application of Lagrange multipliers to verify that the unique element of L 2 
with minimum norm which satisfies the constraints (6.17) must be of the form 

n 

^Ori) = ^ (6.18) 

k=i 

for some {AfcjJLj C €. The constants A*. can be determined b}' substituting (6.18) into 
equations (6.17) and solving the resulting n x n system. The system is of the form MX = d 
where M = is an n by n matrix, A is the n vector (Ai, . . . , A„) T and d is an n vector 
(d(ai), . . . , d(a n )) T . The entries of Ad are given by 


/ oo poo 

C{(x\ )Cj (rq ) dxi = / c(xi - a ?: )c(: r 2 - a 5 ) dx x . (6.19) 

■ oo J — oo 

The matrix M is clearly Hermitian and in fact is always invertible. To see this, suppose we 
can find some vector A^O with MX = 0. Then X T MX = 0 and we conclude that 

/OO ^ 

c ( :r i “ Oi)Xic(xi — cij)Xj dxi = / 

00 i,j ' ~ 

This implies 


X i C ( X l ~ °i) 

i— 1 


dx.\ — 0. 


XjC(Xi - CLj) = 0 . 

3 = 1 

Fourier transform both sides and use the basic properties of the Fourier transform to obtain 


c(y)^ZXje ia > y = 0 . ( 6 . 20 ) 

j = i 

The functions fi(y) — e taiy are linearly independent for distinct a i: and analytic, so that 
f(y) — i has isolated zeroes. Based on equation (6.20) we conclude that c(y) = 0 

in Z/-(1R), contradicting Lemma 4.2. Therefore M must be invertible. This inversion proce- 
dure thus always produces a unique estimate of S if the measurement locations are distinct. 


We can also ‘'solve” the inverse problem by choosing the unique function S which satisfies 
equations (6.1 ( ) and has minimal norm in a weighted T 2 space L 2 (IR,) with norm defined by 
the inner product 

roo l 

< f,9 > 5 = J jixMx^dX! 
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where S(x j) is some real- valued non-negative function on 1R. In this case, we have 


/ °° TOO 1 

S(xi)ci(xi) dxi = / 5(xi)c t -(.ri)5(ari )— — - dx x 

-oo 7— oo o(aq) 

where we must assume that S = 0 wherever <5 = 0. Thus the integral is understood to be 
taken only over that set where 6 is non-zero. Equations (6.17) now take the form 

<S,CiS> s =di ' (6.21) 

and the minimal norm solution is of the form 

n 

S(.Ti) = S(x i) J2 X i c i( x i)- (6.22) 

i~l 

The idea is to choose <S(aq) to have the same general form as S'(aq), and so incorporate a 
priori information into the reconstruction based on (6.22) by forcing it to have the same 
general form. For example, if we know that S is supported in the interval [—6, b\ w T e can 
choose S(x) = 1 on [-—6, ft] and S(x) = 0 elsewhere. The optimal estimate of S becomes 

n 

S(x) = Y[-6,6] X i c i( x ) 

i = 1 

where X[-6,6] is the characteristic function of the interval [—6, b] and where the A, are found 
bv solving 

it, i^J_ b c i( x )cj( x ) dxj = dj 

for j = 1 to n. 

Whether we use a priori information or choose a uniform weighting on L 2 , the recon- 
struction procedure is as follows: Compute the matrix M defined by equation (6.19) and 
■‘measure ’ the data d { at the corresponding points aq = on the top surface. Solve the 
system MX — d to obtain Ai, . . . , A„ and then compute an estimate of S(aq) from equation 
(6.18) or (6.22). The stability of the finite data inversion is thus determined by the nature 
of the matrix M , and specifically, of its inverse. We can quantify the stability of the finite 
data inverse problem by studying the conditioning of the inverse of M in various situations. 
This is done in the following section by studying the singular values of the matrix M. 
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7 Numerical Experiments 


We will now examine the finite data version of the inverse problem bj' using the previouslv 
described inversion procedure. In this section we apply the procedure to simulated data 
sets, both with and without noise. Our main focus is to examine the stability and resolution 
of back surface estimates with respect to various experimental parameters, specifically, the 
frequency of the input heat flux and the distribution of the measurement locations along 
the top surface of the sample. We also demonstrate how a priori assumptions about the 
nature of the corrosion can be incorporated into the inversion procedure, and the effects 
such assumptions have on stability and resolution. 

There are a few points worth mentioning before we present the numerical results. The 
estimates of 5 constructed using equation (6.18) lie in H°°( IR) and so the graph x 2 = S(x i) 
would make sense as a curve in IR 2 , except that the estimate will usually be complex-valued. 
This is not surprising, for a complex-valued 5 makes perfect sense in the linearized version 
of the problem on wdiich the inversion procedure is based. Estimates of 5, especially in the 
presence of noise or the linearization error, will almost certainly have non-zero imaginary 
part. However, for those S of “small” norm the linearized problem accurately reflects the 
full non-linear (with respect to 5) direct problem and so the estimate of 5 should be “mostly 
real” , that is, it should have a relatively small imaginary component. This is indeed the case. 
A physically meaningful estimate of the true back surface x 2 — S^) can be provided by 
either dropping the imaginary component or taking the modulus of the estimate. We choose 
the latter. 

In the examples that follow we generate simulated test data using the full direct problem 
(2.2) with heating g(x) = 1. The direct problem is solved by converting it into a boundary 
integral equation wdiich is then solved numerically. The boundary integral formulation leads 
to a second kind Fredholm equation (explained below). The fact, that g is not compactly 
supported nor even L 2 presents a minor problem. This can be fixed by simply subtracting 
off the function u () satisfying the direct problem with 5 = 0. Recall that u 0 can be found 
explicitly by solving equation (3.5). The function v = u - u 0 then satisfies the boundary 
value problem 
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Au — iuov = 
dv 

d^ 

dv 

dv 


0 in 


0 on x 2 = 1, 


duo 

dv 


on x 2 = 5(rri). 


If 5'(;ri) is sufficiently rapidly decreasing then the Neumann data in the above boundary 
value problem is L 2 and the integral equation formulation becomes 


1 , , f dG 

-v{x) + / — 

2 JdQo dv y 


(x, y)v{y) da y 



dcr y 


where da y is surface measure on dfl 0 and G(x,y ) = V(\x - y\). The integral equation is 
solved using Nystrom’s method (see [2]) with appropriate quadrature rules on the interval 
(— 00 , 00 ) that arises when the boundary of fl 0 is parameterized. The boundary integral 
approach is exceptionally fast and accurate and provides estimates of the solution only on 
the bound an' of S2 0 > the only place that we need the solution. 

To illustrate the general procedure and to show that the inversion algorithm provides 
reasonable estimates, we begin with a simple example. We apply the inversion procedure to 
data generated using the back surface 


S(x) = 


e -(x-l ) 2 


10 


+ 


e -tr+2) 2 /2 

5 “ 


We use a heating frequency of u — 1. As a first step the function c(x ) defined in Lemma 4.1 
is computed by solving the system (4.10) with <2 = 0 numerically and then computing 
^0 — “2r + t j o, again via Nystrom’s method. These quantities are precomputed and stored 
rather than re-computed every time they are needed, for the function c(x) depends only on 
^ 0 ’ not S. The same is true of the matrix M . The entries of the matrix are computed 
by using equation (6.19) and an appropriate quadrature rule 011 (— 00 , 00 ). I 11 the example 
below the data is computed using the full direct problem. The temperature data vector d 
is computed at 21 equally spaced points on the top surface, x 1 = a { where a { = — 5 + | for 
7 — 0 1° 20. We then invert the 21 x 21 system MX = d to find A and return an estimate of 
S via equation (6.18). The estimate of S is computed at a suitable number of points on the 
lange of interest, in this case from — 5 to 5. The reconstruction is shown in Figure 2. The 
dotted line is the actual function S(x) and the solid line is the reconstructed version. 
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Stability 

Of particular interest is the sensitivity of the inversion procedure with respect to various 
experimental parameters, e.g., heating frequency and measurement locations. The first task 
is to quantify the stability or conditioning of the finite data inverse problem. One sensible 
way to do this is to perform a singular value decomposition on the matrix M defined by 
equation (6.19) and examine the magnitude of the singular values. When the singular values 
are small the inversion of MX = d magnifies small perturbations in d. Put another way, 
small singular values mean that relatively large changes in S (and so in A) produce relatively 
small changes in the data, so that perturbations in the back surface are “hard to see.” Our 
goal in choosing experimental parameters is therefore to make the singular values of M as 
large as possible, within certain limits. We should remark that one can define the condition 
number of M as the ratio of the largest to smallest singular values and then attempt to 
quantify stability using this single number. This is not always good approach in the present 
setting, as later examples will show. 

Let us begin by examining how the stability of the inversion procedure depends on the 
locations of the temperature measurements on the top surface. In the following examples 
we fix the heating frequency at uj — 1 and take measurements of the resulting temperature 
at 21 equally spaced locations on the interval [—a, a] for several values of a. The resulting 
measurement locations are therefore of the form a* = — a + ^ for * = 0, . . . , 20. In each case 
the matrix M is computed and a singular value decomposition is performed. Let the 
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— — — a = 10.0 

a = 5.0 

a = 3.0 

a = 2.0 

a = 1 .0 


Figure 3: log 10 |a*-| versus i for various values of a. 


singular values of M be denoted by a,-, i — 1 to 21, arranged in descending order. In Figure 3 
we plot the quantity log 10 |ct* | versus i for the cases a = 1,2, 3, 5, 10. 

It is apparent that as the measurement locations become more spread out (as a gets 
larger) the singular values decay more slowly and hence the inversion procedure becomes 
more stable. In light of Theorem 5.2 this is not surprising. When the measurement locations 
are close together we are able to resolve higher spatial frequencies in the data and so we 
are able to estimate higher frequencies in the Fourier decomposition of 5. But according to 
Theorem 5.2, these are exactly the portions of 5 that are difficult to reconstruct — they are 
heavily damped out in the data. The finite data version of the problem reflects this, with a 
full 6 orders of magnitude variation for the smallest singular values between the cases a = 1 
and a = 10. 

Another way to look at the stability of the various experimental configurations is to sup- 
pose that we have an "error magnification tolerance” E, and that in the inversion procedure 
we disregard all singular vectors whose singular values are less than This idea has been 
used in studying the stability for the impedance imaging problem [8], The inversion proce- 
dure is then stabilized at the expense of rendering those components of S lying in the span 
of the corresponding functions invisible. Figure 4 shows the number of singular values of M 
which satisfy > ~ versus log 10 (F T ) for E from 1 to 10~ 9 . As in the previous examples, the 

matrix M is 21 x 21 and we use measurement locations on the top surface a, = —a 4- — 

r 1 10 ’ 

i = 0, . . . , 20 for a = 1, 2, 3, 5, 10. The heating frequency is u) — 1. 
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Figure 4: Number of singular values with ^ versus log 10 (j£) for various values of a. 


Figure 4 also makes clear that, as the measurement locations become spread out more 
singular values satisfy The inversion procedure then admits more basis functions, 

presumably improving the fidelity of the reconstruction. In the two cases below we perform 
the actual reconstruction with E = 100 (so only singular values greater than 0.01 are ad- 
missible) and add a small amount of random noise to the data (equal to 10 percent of the 
maximum signal strength). We then perform a reconstruction which omits all basis vectors 
whose corresponding singular values are less than Figure 5 illustrates the case in which 
the measurements locations are equally spaced from — 5 to 5; there are 9 admissible singular 
values. 



Figure 5: Reconstruction of S(x) for 21 measurements on [—5,5], tolerance E = 10 2 


In Figure 6 we take the 21 measurements on the smaller interval [—1,1], which yields only 3 
admissible singular values. 
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Figure 6: Reconstruction of S(x) for 21 measurements on [-1, 1], tolerance E = 10 2 . 

The reconstruction in Figure 6 is noticeably inferior to that of Figure 5, but we have only 3 
admissible basis functions with which to construct S(x). Increasing the value of E to admit 
more basis functions is not successful. Figure 7 illustrates what happens if we take E = 10 4 
with measurements on [—1, 1]. Now 5 singular values are admissible, but the reconstruction 
is overwhelmed by noise. 



Figure 7: Reconstruction of S(x) for 21 measurements on [-1,1], tolerance E = 10 4 . 

The moral seems clear, for maximum stability with a fixed number of measurement 
locations, we should spread the measurements over as large a region as possible. There are 
limits to this approach, however. If we spread out the measurements we do gain stability, 
but we will no longei be able to estimate high frequencies in the Fourier decomposition of 
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S. This is illustrated by Figure 8, where we take 21 noise-free measurements on the interval 
[—10, 10] and estimate S with error tolerance E = 10 2 . In this case all of the singular values 
are admissible. 



Figure 8: Reconstruction of S(x) for 21 measurements on [—10, 10], tolerance E — 10 2 . 

Despite the fact that the inversion is quite stable, our inability to resolve high frequencies 
results in a loss of resolution of small-scale detail in the reconstruction. With regard to the 
distribution of the measurement locations, the reconstruction process involves a compromise 
between stability and resolution of small-scale features. 

In the next series of examples we examine the dependence of the stability on the 
frequency of the applied heat flux. We consider the cases uo — 0.01,0.1,1.0,10.0,100.0. 
In each case we take 21 equally spaced temperature measurements on the interval [—5,5]. 
Figure 9 shows log 10 |cq| versus i for each case. 

w = o.oi 
w = o.i 
w = 1 .0 
w = 10.0 
w = 100.0 


Figure 9: log 10 |a,| versus i for various values of co. 

The figure illustrates that higher frequencies give rise to much smaller singular values. 
As before, it is instructive to consider the case in which we have an error tolerance E and in 
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the inversion process we omit those singular vectors whose singular values are less than 
Figure 10 shows the number of singular values which exceed ~ for E from 10~ 8 to 10 15 , for 
each of the frequencies. 



“5 0 5 10 15 


Figure 10: Number of singular values with > T versus log 10 (£') for various values of to. 

Figure 10 clearly illustrates the situation. For E = 10~ 2 as before, uj = 10 and uj = 100 
have no admissible singular values at all. The reconstruction (if carried out) is identically 
zero. The smaller singular values at higher frequencies are due to the fact that at higher 
frequencies the periodic heating penetrates very little into the sample and becomes more of 
a “skin effect.” As a result very little energy reaches the back surface and even less returns 
to be measured on the top surface; at high frequencies the map taking S into d is essentially 
multiplication by zero. It is interesting to note that while higher frequencies produce smaller 
singular values, the condition number of M for u — 100 is only 10.3, while the condition 
number for uj ~ 0.01 is 706.1. Clearly, though, it ? s not enough to make the condition number 
small. The singular values themselves must be large enough for the inversion procedure to 
be stable in the presence of a fixed noise level. 

While lower frequencies make the inversion process more stable, there are limits to how 
small we can make uj and maintain resolution. Figure 11 illustrates a reconstruction based 
on uj = 0.1. The parameters are otherwise identical to those that were used to produce 
Figtue 5. All of the singular values are admissible. In fact, the smallest singular value is 
0.063. As with the case in which the measurement locations were spread out over [—10, 10], 
we lose resolution at low temporal frequencies for the input flux. 
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Figure 11: Reconstruction of S(x) for 21 measurements on [—5,5], to = 0.1, tolerance 

E = 10 2 . 


Incorporating A Priori Information 

The preceding examples illustrate that the inversion procedure involves a compromise 
between stability and resolution. If the data points are too closely spaced, the inversion 
procedure is unstable. If the data points are too spread out, the inversion procedure becomes 
stable, but resolution is lost; measurements taken far from the support of the defect contain 
little information, because the heat diffuses very rapidly. Variations in the input heating 
frequency give rise to a similar phenomena. How shall we find the “best” experimental 
parameters? One useful possibility is to incorporate a priori information or assumptions into 
the inversion procedure. We will illustrate the idea by examining the problem under the 
assumption that the defect or function S is supported in a known interval. 

In the following examples we assume that the defect being imaged is supported in the 
interval [—2,2]. The only modification to the inversion procedure is that the matrix M is 
computed in accordance with equation (6.21) and the function 5 is estimated using equation 
(6.22). We will study the stability of the inversion procedure with respect to the distribution 
of the measurement locations on the top surface. 

As in the previous cases, we choose measurement locations at x\ — a { on the sample top 
surface, where a z - = —a + -~a for i = 0 to 20. The heating frequency in all cases that follow 
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a = 0.5 
a = 1.0 
a = 2.0 
a = 5.0 
a = 10.0 


Figure 12: Singular values a * versus t for various values of a. 

is a; = 1. Let us begin by examining the singular values of the inversion matrix M for a few 
choices of a. In Figure 12 we plot the quantity log 10 |cq| versus i for a = 0.5, 1.0, 2.0, 5.0, 10.0. 

The figure shows that the best conditioning for the inverse problem occurs at a = 2, 
when the measurement locations are distributed approximately in the same interval in which 
the defect is assumed to be supported. As before, closely spaced locations give rise to 
an ill-conditioned problem. However unlike the previous cases widely spaced nodes also 
result in poor conditioning. When M is computed using equation (6.21) those rows of M 
corresponding to measurement locations far from the support of S are very nearly set to zero 
since the function c(x — a*-) is rapidly decreasing away from a*. 

If an error magnification tolerance E is specified, we can plot the number of allowable 
singular values cv*- > versus 'log 10 (2?) for the different node spacings: 

a = 0.5 

— — — a = 1.0 

a = 2.0 

a = 5.0 

a =10.0 


2.5 5 7.5 10 12.5 15 17.5 

Figure 13: Number of singular values with cv 2 - > Jr versus log 10 (F^) for various values of a. 

As expected, a = 2.0 allows more singular values for a fixed value of E than any other 
choice for measurement spacing. It is useful to look at a few reconstructions based on this 
strategy. In the two cases below we take E = 300 (so only singular values greater than 
are admissible) and add a small amount of random noise to the data (equal to 1 percent of 
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the maximum signal strength). We then perform a reconstruction which omits all singular 
values less than jjj. The function defining the back surface is S(x) = T- e -2(>-fi) 2 + ±-e~ 3(x ~ 1)2 . 
Figure 14 illustrates the first case using a = 2, the best choice according to Figure 13. In 
this case 7 singular values are admissible. 



Figure 14 : Reconstruction of S(x) for 21 measurements on [—2,2], tolerance E = 300. 
For a — 10 we have 4 admissible singular values and the reconstruction shown in Figure 15. 



Figure 15: Reconstruction of S(x) for 21 measurements on [—10, 10], tolerance E = 300. 

The case a — 0.5 also yields 4 admissible singular values and the reconstruction shown in 
Figure 16. 
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Figure 16: Reconstruction of S(x) for 21 measurements on [-0.5, 0.5], tolerance E = 300. 

The actual reconstructions confirm that a = 2 yields the most desirable results. Choosing 
a significantly smaller or larger than the support of S results in decreased stability and/or 
accuracy for the reconstruction. 

Of course, the assumption that S is supported in a given interval should be detrimental 
to the leconstruction if that assumption turns out to be false. In the following case we let 
= T e - 2 (.r+ip + I e -3(*- 4 ) 2 and p er f orm the reconstruction under the assumption that 
5 is supported in the interval [—2,2], We take measurements at 21 equally spaced location 
between —2 and 2, the best case from above, and use an error tolerance E — 300. The result 



Figure 17: Reconstruction of S(x) for 21 measurements on [—2,2], tolerance E — 300. 
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The incorrect assumption obviously introduces errors into the reconstruction, although that 
portion of S which is non-zero in the interval [—2, 2] is still recovered with reasonable accu- 
racy. 

Concluding Remarks. 

In this paper we have investigated the inverse problem of recovering. an unknown bound- 
aiy portion of some object by applying a heat flux to an accessible portion of the boundary 
and measuring the resulting temperature response. We have considered a linearized version 
of the problem and found that the continuous version of the inverse problem, in which one 
has data at every point on the accessible portion of the surface, is extremely ill-posed. In- 
deed, the linearized version requires one to solve a first kind convolution integral equation 
for the unknown surface. The convolution kernel has a Fourier transform which dies rapidly 
at. infinity, and so the inversion is extremely sensitive to the data at high spatial frequencies. 
We performed a variety of numerical studies which show that the ill-posedness is directly 
reflected in the finite data version of the problem, by the rapid decay of the singular values 
of the matrix which governs the inversion process. This ill-posedness depends on a num- 
ber of factors; in particular, the locations of the measurements have a large effect on the 
conditioning of the inverse problem, and these effects mirror the behavior of the continuous 
version. We have also considered the effect of including a priori assumptions in the finite 
data inversion procedure, b}^ weighting appropriate Hilbert spaces in which the solution S 
resides. The inclusion of this information can help in determining the optimal locations for 
measurements on the top surface. 

There are a number of interesting directions we could take from here. In our studies 
we used only the input flux whose magnitude is identically one on the top surface. Similar 
results can be obtained for more general fluxes, and this would allow one to study the effect 
that the input heat flux has on sensitivity and resolution. The fully time-dependent case 
would also be of interest. The procedure presented in this paper would also work for a full 
three- dimensional problem, although qualitatively the results should be the same — the high 
spatial frequencies in the back surface should be difficult to see. 

As mentioned earlier, the inversion process which chooses that function with minimal L 2 
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norm which is consistent with the measured data seems to act like a form of regularization for 
the inverse problem. It would be interesting to examine in what sense this is true, and how 
it relates to more traditional forms of regularization. It is also possible (and not difficult) to 
carry out the same minimization process in higher Sobolev spaces, e.g., HK and thus put a 
higher penalty on functions with oscillations. This too would make an interesting study. 
We would also like to examine conditions under which our inversion procedure is guaranteed 
to converge to the solution of the linearized inverse problem. 
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